home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YDEGEN.C < prev    next >
C/C++ Source or Header  |  1988-11-28  |  4KB  |  154 lines
  1. /********************** YDblRtr.C   for the double rotor *********************/
  2. /*********************** (C) 1986,7,8 by JAMES A. YORKE **********************/
  3.  
  4.  
  5. #include "yinclud.h"
  6.  
  7. double  exp();
  8.  
  9. static double   u11,
  10.                 u12,
  11.                 u21,
  12.                 u22,
  13.                 s1,
  14.                 s2;
  15.  
  16.  
  17.  
  18. /* for double rotor: */
  19.  
  20. static double   Matrix11,
  21.                 Matrix12,
  22.                 Matrix21,
  23.                 Matrix22;
  24. static double   cc1,
  25.                 cc2;
  26.  
  27.  
  28.  
  29. /*********** Degenerate case with hitting time = infinity ***************/
  30. startDegenerateDblRotor() {
  31.  /* calculate parameters for the Degenerate double rotor map */
  32.     double  mass1,
  33.             mass2,
  34.             length1,
  35.             length2,
  36.             damping1,
  37.             damping2;
  38.  /* above parameters mostly to be set from menu */
  39.     double  inertia1,
  40.             inertia2;
  41.     double  discrim,
  42.             unorm;
  43.     mass1 = C1;        /* mass1 */
  44.     mass2 = C2;        /* mass2 */
  45.     length2 = C4;        /* length2 -- length1 is set below */
  46.     damping1 = C5;        /* damping1 & damping2 are damping  */
  47.     damping2 = C6;
  48.  
  49.  /* Let the length of the first rod and the inertias be set to values that
  50.     allow for a symmetric set of ODEs for the system */
  51.     length1 = length2 * sqrt(mass2 / (mass1 + mass2));
  52.     inertia1 = (mass1 + mass2) * length1 * length1;
  53.     inertia2 = mass2 * length2 * length2;
  54.  
  55.  /* compute eigenvalues for revised map */
  56.     discrim = 0.5 * sqrt(damping1 * damping1 + 4.0 * damping2 * damping2);
  57.     s1 = 0.5 * damping1 + damping2 + discrim;
  58.     s2 = 0.5 * damping1 + damping2 - discrim;
  59.     u11 = damping2 - s1;
  60.     u12 = damping2;
  61.     unorm = sqrt(u11 * u11 + u12 * u12);
  62.     u11 = u11 / unorm;
  63.     u12 = u12 / unorm;
  64.     u21 = damping2 - s2;
  65.     u22 = damping2;
  66.     unorm = sqrt(u21 * u21 + u22 * u22);
  67.     u21 = u21 / unorm;
  68.     u22 = u22 / unorm;
  69.  
  70.  /* Compute matrix M(T) for theta dot */
  71.  
  72.  /* Compute matrix M(T) for theta  */
  73.     Matrix11 = -(u11 * u11 * (-1) / s1 + u21 * u21 * (-1) / s2);
  74.     Matrix12 = -(u11 * u12 * (-1) / s1 + u21 * u22 * (-1) / s2);
  75.     Matrix21 = -(u11 * u12 * (-1) / s1 + u21 * u22 * (-1) / s2);
  76.     Matrix22 = -(u12 * u12 * (-1) / s1 + u22 * u22 * (-1) / s2);
  77.  /* miscellaneous constants -- rho is forcing and is only in iterated map */
  78.     cc1 = length1 / inertia1;
  79.     cc2 = length2 / inertia2;
  80. }
  81.  
  82. initDegenerateDblRotor() {
  83.     extern int      startDegenerateDblRotor();
  84.     extern int      DegenerateDblRotor();
  85.  
  86.  
  87.     X_upper = pi;        /* x scale */
  88.     X_lower = -pi;
  89.     Y_upper = pi;        /* y scale */
  90.     Y_lower = -pi;
  91.  
  92.     vec_dim = 2;        /* the dimension of the Lyapunov vectors =
  93.                    phase space dim */
  94.     num_lyap = 0;        /* the number of Lyapunov numbers to be
  95.                    computed <= vec_dim */
  96.     lyapzero = 2;        /* y[lyapzero] is the zeroth coord of the
  97.                    zeroth lyapunov vector */
  98.     dim = lyapzero + num_lyap * vec_dim;
  99.  /* needed for rungekutta */
  100.  
  101.  /*  mass1,2 = C1,C2;     length2 = C4; length1=C4*sqrt(C2/ (C1+C2));
  102.     damping1,2 = C5,C6      hit time=infinity   */
  103.     C1 = 1.;
  104.     C2 = 1.;
  105.     C4 = 1.;
  106.     C5 = 1.;
  107.     C6 = 1.;
  108.  
  109.     rho = 2.5;        /* forcing strength */
  110.  
  111.     init_process = startDegenerateDblRotor;
  112.     map = DegenerateDblRotor;
  113. }
  114.  
  115.  
  116. DegenerateDblRotor() {
  117.     double  moduloAB();
  118.     double  rhosiny0,
  119.             rhosiny1,
  120.             rhocosy0,
  121.             rhocosy1;
  122.     int     i;
  123.     int     base;
  124.     double  y_temp[6];
  125.     rhosiny0 = rho * sin(y[0]);
  126.     rhosiny1 = rho * sin(y[1]);
  127.     y_temp[0] = y[0] + Matrix11 * cc1 * rhosiny0 + Matrix12 * cc2 * rhosiny1;
  128.     y_temp[1] = y[1] + Matrix21 * cc1 * rhosiny0 + Matrix22 * cc2 * rhosiny1;
  129.  
  130.  
  131.     if(num_lyap > 0) {
  132.         rhocosy0 = rho * cos(y[0]);
  133.         rhocosy1 = rho * cos(y[1]);
  134.     }
  135.     for(i = 0; i < num_lyap; i++) {
  136.         base = lyapzero + vec_dim * i;
  137.         y_temp[base]
  138.             = y[base] + Matrix11 * cc1 * rhocosy0 * y[base]
  139.             + Matrix12 * cc2 * rhocosy1 * y[base + 1];
  140.         y_temp[base + 1]
  141.             = y[base + 1] + Matrix21 * cc1 * rhocosy0 * y[base]
  142.             + Matrix22 * cc2 * rhocosy1 * y[base + 1];
  143.     }
  144.     for(i = 0; i < lyapzero + vec_dim * num_lyap; i++)
  145.         y[i] = y_temp[i];
  146.  
  147.  /* now bring y[0] and y[1] into -pi,pi etc. */
  148.     y[0] = moduloAB(y[0], -pi, pi);
  149.     y[1] = moduloAB(y[1], -pi, pi);
  150. }
  151.  
  152.  
  153.  
  154.